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A Maximally Localized Wannier Function Approach 
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We present an ab initio inelastic quantum transport approach based on maximally localized 
Wannier functions. Electronic-structure properties are calculated with density-functional theory in a 
planewave basis, and electron-vibration coupling strengths and vibrational properties are determined 
with density-functional perturbation theory. Vibration-induced inelastic transport properties are 
calculated with non-equilibrium Green's function techniques, which are based on localized orbitals. 
For this purpose we construct maximally localized Wannier functions. Our formalism is applied to 
investigate inelastic transport in a benzene molecular junction connected to mono-atomic carbon 
chains. In this benchmark system the electron-vibration self-energy is calculated either in the 
self-consistent Born approximation or by lowest-order perturbation theory. It is observed that 
upward and downward conductance steps occur, which can be understood using multi-eigenchannel 
scattering theory and symmetry conditions. In a second example where the mono-atomic carbon 
chain electrode is replaced by a (3, 3) carbon nanotube, we focus on the non-equilibrium vibration 
populations driven by the conducting electrons using a semi-classical rate equation. 



I. INTRODUCTION 

Molecular electronic devices have been intensively 
studied for the past decade, having been regarded as can- 
didates to replace silicon-based electronic devices. Since 
Aviram and Ratner proposed the concept of the first 
molecular rectifier in 1974*, a variety of molecular devices 
have been suggested. In addition, thanks to advancement 
in experimental fabrication and measurement techniques, 
electronic currents through molecules have been experi- 
mentally measured and investigated 2 - - — . Despite these 
efforts, there are still many issues in the practical real- 
ization of molecular electronic devices. 

Understanding interactions between conducting elec- 
trons and molecular vibrations is one of the key issues 
to address for future applications. Vibrational excita- 
tions due to the scattering of conducting electrons can 
change molecular configurations and attachment geome- 
tries, affecting the functionality and performance of the 
molecule-based devices or in turn backscatter electrons, 
impeding flow. The worst scenario is that local heating 
effects might ultimately break down the junction. 

Vibration-induced inelastic transport has been theo- 
retically investigated following two directions. The first 
one is to study simple model Hamiltonians, e.g. a single 
electronic level coupled to a single phonon mode, which 
is known as the Anderson-Holstein model- - — . Based on 
this simplified model, many novel and interesting trans- 
port properties have been predicted and investigated. 
However, the model used in this approach is too sim- 
plified to provide detailed and accurate theoretical data 
that quantitatively explain experimental results. 

On the other hand, a quantitative and computational 
approach based on density functional theory (DFT) offers 
the chance to describe a realistic system accurately with- 



out any adjustable parameter s 17 ' 18 . In particular, using 
DFT one can calculate equilibrium geometries, electronic 
couplings, normal modes, and electron- vibration interac- 
tions. 

Transport theories can then be combined with DFT, 
and several approaches have been proposed 1 ^ - —. In 
particular, a non-equilibrium Greens Function (NEGF) 
formulation, in combination with DFT, and commonly 
called DFT-NEGF, has been widely used in ab initio 
quantum transport problems 1 ^ - — . This approach is more 
powerful than other methods in that it can tackle not just 
the emergence of electron- vibration interactions, but any 
other type of interactions. DFT-NEGF has been success- 
fully applied to elastic quantum transport for both zero- 
bia o 22 ' 23 and finite-bias cases 1 ^ - — , and recently it has 
been extended to include interaction effects like electron- 
vibration interaction s 40 ' 42 . 

DFT-NEGF requires to use an atomic-like localized 
basis, because the device should be spatially divided into 
two electrodes and a molecular conductor. For this rea- 
son most of DFT-NEGF calculation packages have been 
implemented using localized basis set. However, from a 
computational viewpoint, it is known that a planewave- 
based DFT calculation can provide an accurate descrip- 
tion of electronic states systematically, and in particular 
can describe electronic states which have considerable 
spread in vacuum, where localized orbitals are absent. 
Furthermore, while basis functions used in localized-basis 
calculations are determined depending on types of atoms 
and the chemistry of the system, a planewave basis can 
describe any given system without making any further 
assumption. 

However, periodic-boundary conditions planewave cal- 
culations are not suitable to DFT-NEGF calcula- 
tions. For this, maximally localized Wannier functions 
(MLWF), as proposed by Marzari and Vanderbilt 31 , 



2 



provide the formal and algorithmic formulation for a 
transformaton between delocalized and localized orbitals. 
Since the Wannier transformation is an exact unitary 
mapping, one can construct a minimal set of atomic-like 
localized functions within an energy window of interest 
without losing the accuray of planewave-based DFT cal- 
culations, and a MLWF approach to quantum transport 
has been very successfully applied to zero-bias quantum 
conductance calculation s 22 ' 23 . The next step to develop a 
MLWF approach to quantum transport 2 ^— is to include 
interaction effects on transport properties. In this paper 
we focus on extending MLWF-based ab initio quantum 
transport calculations to investigate electron-vibration 
interaction effects on molecular junctions. 

This paper is organized as follows. In Sec Hi] we re- 
view first-principles electronic structure calculations, es- 
pecially focusing on (1) density- functional perturbation 
theory to calculate vibrational properties and electron- 
vibrational interaction, and (2) the transformation of 
electron- vibration interactions from a plane- wave basis to 
a maximally localized Wannier function basis. In Sec lIIII 
a quantum transport theory based on non-equilibrium 
Green's functions and diagrammatic perturbation the- 
ory is presented. First we discuss electrode-conductor- 
electrode system partitioning and the calculation of lead 
self-energies. Then we review diagrammatic expansion 
schemes for electron-vibration interactions in the self- 
consistent Born approximation and lowest-order pertur- 
bation theory. Inelastic transport properties such as 
finite-bias electronic current and power transfer are cal- 
culated within the Meir- Wingreen transport formalism 39 . 
Finally we present a semi-classical rate equation to de- 
termine nonequilibrium vibrational populations in the 
presence of interactions with conducting electrons and 
coupling to bulk vibrations in the electrodes. In SecHVl 
application results and further analysis are presented. 



II. ELECTRONIC STRUCTURE METHODS 

A. Vibrational Properties: Density-Functional 
Perturbation Theory 

The electronic structure of a given system is calculated 
within the DFT framework. The ground state charge 
density and Bloch wavefunctions are determined by solv- 
ing the Kohn-Sham equations 2 ^, 
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where 



V KS 



n(r') 



-dv' + v xc (r) + V lon {r). (2) 



Here v xc = H^s. is the exchange-correlation potential, 
and Vi on (r) is the ionic core potential described by pseu- 
dopotentials. We solve the Kohn-Sham equations in 



a planewave basis, as implemented in the QUANTUM- 
ESPRESSO distributional. 

Vibrational properties and electron- vibration interac- 
tions are determined within density-functional perturba- 
tion Theory (DFPT) 2 -. Vibrational spectra and the cor- 
responding normal modes are obtained from the inter- 
atomic force constants, which are the second derivatives 
of the Born-Oppenheimer total energy surface i?({R}) 
with respect to displacements of the ions, 
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where i(J) represents the i(j)th atom of the unit cell, and 
a(/3) indicates the Cartesian component of the displace- 
ments. While C^°"sj * s ^ ne contribution of the ion-ion 
interaction potential, C^ e ^ is the second derivative of 
electron-electron and electron-ion interactions. From the 
Hellmann-Feynman theorem, one can obtain 

n(r)- 
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The key quantity to obtain interatomic force constants 
is the linear variation of the charge density n(r) with re- 
spect to ionic displacements. In the Kohn-Sham map- 
ping, the linear response of the charge density Ari(r) 
can be calculated from the Kohn-Sham orbital variation 
\A?pi). From first-order perturbation theory, one can de- 
rive the equation for \Aipi)2L: 
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where 
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dv' 



dv r 



dn 



An(r). 
(6) 

Here Aei = (ipi\AVKs\i>i) is the first-order correction to 
the Kohn-Sham energy eigenvalue £j. Note that these 
equations should be self-consistently solved since the 
equations that determine \Aipi) depends on the linear 
response in the charge density An(r). The vibrational 
frequencies uj of the system can be calculated by solving 
the following eigenvalue equation: 
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The electron-vibration interaction can be written in a 
second quantized form as follows^: 

n e l- mb = Sk+£k °k+q a k (&LqA + K>) > ( 8 ) 



kqA ran 



where a^™ (a£) is the electron creation (annhilation) op- 
erator for Bloch state \ip n u)- Similarly 6 qA (bq\) is the 
creation (annhilation) operator for the phonon of the vi- 
brational mode A with the energy w q A at wave vector 
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q, and <?k_|_' q k * s ^ ne electron-phonon coupling matrix el- 
ement. Once the linear response of the charge density 
is computed from DFPT, the electron-phonon coupling 
matrix can be calculated from the derivative of the self- 
consistent Kohn-Sham potential Avks as follows: 



qA,mn 

■9k+q,k — 




where ip^ n is the nth Kohn-Sham orbital wavefunction at 
wavevector k. is the response of the self-consistent 

Kohn-Sham potential with respect to the phonon mode 
A at wavevector q. 

Because we use a large supercell that contains a con- 
ducting molecule and two electrodes, T-point sampling 
can be used safely. Thus, by dropping wavevector indices 
for electrons and vibrations, the interaction Hamiltonian 
can be written in the simpler form: 

H el . mb = E E 9x n ala n (b{ + b x ) , (10) 

A ran 

where 

9x n =(^A ' <VUAi4 s |Vn>. (11) 



B. Maximally Localized Wannier Functions 

As discussed in the introduction, DFT calculations 
based on a planewave basis can provide a very accurate 
description on the electronic structure of the system, in 
particular in comparison with a localized basis set. How- 
ever, since Bloch orbitals are intrinsically delocalized, 
they are not suitable to quantum transport calculation 
based on Green's functions, in which spatial seperation 
between electrodes and the conductor is required in the 
Hamiltonian description. In this work maximally local- 
ized Wannier functions (MLWFs ) 31 i 32 are used in order to 
transform Bloch wavefunctions into localized functions. 
When there are N isolated bands, WFs can be defined 
as 

N 

k rn 

where R is the Bravais lattice vector and the Wannier 
rotation matrix is unitary. In the MLWF construc- 
tion, f/( k ) is determined by minimizing the mean square 
spread of the resulting Wannier functions, defined as 

N 

n 
N 

= ^2 [(w„o|r 2 |w n0 ) - (w„o|r|wno) 2 ] . (13) 
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Note that by construction MLWFs are strongly local- 
ized in real space. Since the Wannier transformation 
is an exact unitary mapping, the numerical accuracy of 
planewave-based DFT calculations is preserved in the 
transformation. 

For a metallic system where all energy bands are con- 
nected together, one needs to perform a disentangle- 
ment^ in which a maximally-connected ./V-dimensional 
subspace is extracted from the entire entangled manifold: 

ic p k)= E (1 4 ) 

(k) 

here, N^ n is the number of entangled bands in a desired 
energy window at a wavevector k. For details, we refer 
to Ref^ 2 .. 

Using Eq. (fT2)) and Eq. (fT4]) . one can readily obtain the 
electronic Hamiltonian and electron- vibration interaction 
Hamiltonian in the MLWF basis. In particular, we pro- 
vide an explicit transformation of Hamiltonian from the 
Bloch representation to Wannier one for the T-point sam- 
pling case, which is used in large supercell calculations. 
Under the Wannier transformation, 

K> = E u in uft\i>j) = E ( udisu ) jn l^>> ( 15 ) 

i,3 3 

the electronic Hamiltonian and electron- vibration inter- 
actions read 

7~Le ^ y T^mnCrn^n (1^) 

Uel- mb = E E M mn4c« (b{ + & A ) , (17) 

A m,n 

where 

W mn = (u> m \H e \u n ) (18) 

= J2(u dls u)* m mne\^)(u d -u) jn (19) 

ij 

M^ n = ^-(cj m \AV x \u; n ) (20) 

ij 

Here cj„(c m ) are the electron creation (annihilation) op- 
erators for Wannier function \u) m ). 

III. QUANTUM TRANSPORT THEORY 

A. System Division 

A two-terminal system can be divided into three parts: 
the left electrode (L), the central region (C), and the 
right electrode (R). The central region is defined as the 
region that contains the molecules, defects, or nanoscale 
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conductors of interest. With a localized orbitals basis the 
corresponding electronic Hamiltonian has the following 
structure, 



H L U LC 
Ue = | Ucl U c H C R 
H RC Ur 



where Ucl = U\ c and Ucr = U\ 



rc- 



U 



L(R)C 



(22) 



is the 



coupling matrix between the left (right) electrode and 
the central region. Note that the central region does 
not necessarily include only the conducting molecule, 
but, as noted in Eq. (1221) . should be large enough to 
make the direct coupling between the electrodes zero, 
Ulr — Url = 0. This can be achieved by including 
some surface atoms of the electrodes into the central re- 
gion and defining them as an extended molecule. For 
the elastic quantum conductance calculation this condi- 
tion is sufficient to select the device region. However, 
in DFT-NEGF it is assumed that the conducting elec- 
trons are scattered by the molecular vibrations only in 
the device region, while the electrons in the electrodes 
are treated as noninteracting quasi-particles. Therefore 
the central region should be large enough to guarantee 
that the electron-vibration interactions converge to zero 
outside it. 

Once the electronic Hamiltonian U e has been obtained, 
the non-interacting retarded Green's function Gq(s) for 
the conductor region can be derived as follows: 



where 



Z r L = H CL - 



l 



-H LC 



(23) 



(24) 



(e + irj) Il-Hl' 
Vr = n C R , - . * v-Ubc- (25) 

(£ + Vn)X R - Ur 

Here, S^/i? are * ne ^ ea< ^ self- energies, which describe the 
effects of the coupling between the conductor and the 
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FIG. 1. Two-terminal geometry used in quantum transport 
calculation. The device region consists of a vibrating molecule 
and some of the bulk electrode relaxed surface layers. 



leads, r\ is an infinitesimal positive number, and X the 
identity matrix. In case of a non-orthogonal atomic ba- 
sis set calculation, X should be replaced by the overlap 
matrix. However, since MLWFs are orthonormal by con- 
struction, we can safely use the identity. Up to now 
electron-vibration interactions have not been included in 
the Green's function. 

Since Ul/r ar e semi-infinite it would be impossible 

to calculate directly [(e + irj)X L / R — %l/r\ ■ How- 
ever, we can use the fact that the localized orbitals of 
the central region have nonzero overlap only in a fi- 
nite range close to the conductor. If Vu R \ denotes 
a projection operator onto the subspace of the left 
(right) electrode that has finite coupling with the cen- 
tral region, H LC ( R c) = 'Pl(r)U LC (rc) and Kcl(cr) = 
H-cl{CR)Pl(r)- Therefore, what we need to calculate is 

'Pl(R) [(e + if)) X L ( R ) - U L (R)] 1 "Pl(r) , which is the sur- 
face Green's function. This can be obtained by the effi- 
cient and fast iterative method first suggested by Sancho 
et al^-36. 



B. Interacting Green's Functions 

The full Green's function with the electron-vibration 
coupling can be obtained from diagrammatic perturba- 
tion theory based on non-equilibrium Green's functions, 
known as the Keldysh formalism. From the theory one 
can obtain four kind of interacting Green's functions: 
G r ' a (e) and G>(e). They satisfy the Dyson equation and 
Keldysh equation respectively, 



Qr/a _ Q r / a _|_ Q r / a Yf/ a Q r / a 



G$ = G r (>:,; • >:',; • \ .:, ) 



(26) 

The lead lesser and greater self-energies and are 
£<(e) = ? r Q ( £ )/ Q ( £ ) (28) 
E>( E ) = *T Q (e)(/ a (e)-l), (29) 

where T a = i (£> — £<) and f a is the Fermi-Dirac dis- 
tribution with chemical potential /i Q . In principle the 
electron-vibration self-energy £^{j and S^, contain all 
possible diagrams that satisfy Feynman rules. In this 
work the electron-vibration self-energies £l{f and sj if) 
are calculated from the lowest order diagrams in the di- 
agrammatic expansion. The lesser and greater electron- 
vibration self-energies Y>> ib (e) are 



/OO J I 
(e-e')M x G$ {e')M{m 

where DJ~(e) are the zeroth-order unperturbed correla- 
tion functions for the Ath vibrational mode, 



D>(e) = -2m [N x 5 (e T u\) + (N x + 1) S (e ± w A )] . 



(31) 
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Here N\ is the vibronic occupation number. From the 
useful relation between the self-energies S r — S a = S > — 
T, < and the Kramers-Kronig relation one can readily ob- 
tain 
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'-H[X> ib -Z< b ], (32) 



where £[[■■■] is the Hilbcrt transformation. Equations 
([50)1 and ((321) together with Eq. (gSJ) and (JITJ) constitute 
a set of coupled equations. These equations can be solved 
in the self-consistent Born approximation (SCBA) 4 ° i 42 ' 43 . 
However, the SCBA is computationally very demand- 
ing, and in case of weak electron-vibration coupling, the 
lowest-order perturbation theory (LOPT) may be a good 
approximation. In this approximation, the full corre- 
lation function G^ in Eq. ([30|) is replaced by the non- 
interacting Gq': 
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Within the same order of electron-vibration coupling 
A4\, the Green's functions are calculated as follows: 



G r/a « Go /a + C 
- G r /a + C 
G$ « GS (Ef 
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where G^ = G 



r/a _ f~,r/a v r/a,(2)gr/a 



''vib 



C. Current Calculations 

Once the full interacting Green's functions and self- 
energies are calculated, physical quantities of interest can 
be obtained. In order to calculate transport properties, 
we adopt the Meir-Wingreen transport formula, which 
is widely applied to mesosopic and nanoscale transport 
problems 3 ^. The net electric current I a entering the elec- 
trode is 
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f°° de 
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J-oo 27T 



/ -Tr[E>(e)G<( £ )-E<( £ )G>( £ )](37) 

J-oo L ~ 



where N a is the electron number operator of the a elec- 
trode. The correlation function G^ in Eq. (l3"Tl) is calcu- 
lated in the SCBA and the LOPT, as introduced in the 



previous section. It should be noted that these approxi- 
mations for the Green's functions satisfy the current con- 
servation condition, J2 a lot = 0- ^ ^ s tempting to solve 
Eq.® and together with Eq.® and Ea.QSty, an 
approach known as the first Born approximation (1BA). 
Although this approximation is much easier to calculate 
than the SCBA, the 1BA does not generally guarantee 
current conservation. 

Similarly, one can calculate the net energy current col- 
lected in the electrode, 

Pa = \f |^Tr [£>( £ )G<( £ ) - £<( £ )G>( £ )] • 

(38) 

According to energy conservation, the total sum of the 
energy current to the electrodes and the local vibrations 
should be zero, i.e. Pl + Pr + Pvib = 0. From this 
condition, the net power transferred to local molecular 
vibrations is given by 



A 

= E 



de 
2^ 



£ Tr 



Z> biX (s)G<(s)-E< biX (s)G>m 



D. Vibrational Population 

Due to electron-vibration interactions, conducting 
electrons can exchange their energy with molecular vibra- 
tions by absorbing or emitting local vibrational quanta. 
Therefore the vibrational population N\ in Eq. ([31]) does 
not in general follow a Bose-Einstein distribution. In or- 
der to describe a non-equilibrium vibrational population, 
a semiclassical rate equation has been suggeste d 42 ! 44 , 



dt 



Px 



7a (Na - n B {tuox)) ■ 



(40) 



The first term on the right hand side represents the net 
molecular vibrations excited by the conducting electrons. 
The second term describes a heat dissipation process that 
makes local vibrations equilibrated to the Bose-Einstein 
distribution at the electrode temperature. Molecular vi- 
brations are not isolated, but are mechanically coupled 
to the surrounding (e.g. the bulk phonons of the elec- 
trodes). Due to this coupling, local vibrations in the 
molecule decay into bulk phonons, and 7a represents this 
decay rate. Note that here we restrict ourselves to con- 
sider only the harmonic coupling to bulk phonons. Es- 
sentially the harmonic coupling describes one local vibra- 
tion to one phonon transitions, which mean that a local 
vibration decays by exciting one phonon mode in the elec- 
trodes. Following Ref.(— ), the first term in Eq. (|40|) can 
be decomposed into absorption and emission processes, 



dt 



Nx = {Nx + l)E x - N X A X - 7A {N x - n B (}iu x )) . 

(41) 
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Here the absorption and emission rates are expressed as 
follows: 

9 r°° rip 

A ^ = T\ ? Tr [ M XG< ( £ ~ ^x)M x G> (e)] (42) 

The steady state solution for the vibrational populations 
can be immediately obtained from Eq. (|4"Tj) . 

-4a + 7a - E x 

Note that this steady state solution exists only when 
A\ + 7a > -Ea- Otherwise, Eq.(j4"Tj) gives an exponen- 
tially increasing population as a function of time, which 
implies a vibrational instability. 

In order to calculate vibrational decay rates, let us 
consider a generic case where localized vibrations are 
harmonically coupled to bulk phonons (heat reservoir). 
Localized molecular vibrations and bulk phonons are de- 
scribed by collections of harmonic oscillators in the sec- 
ond quantized form: 

JT S = X>A + 5) (45) 

H B =J2tek(b{b k + ^j, (46) 

where b\(b\) and b~k(b\) are annihilation (creation) op- 
erators for localized vibrations and bulk phonons respec- 
tively. When the Ath molecular vibrational state is occu- 
pied, the decay rate at which the vibration disappears by 
exciting bulk phonons can be calculated by using Fermi 
golden rule^ 5 - - — , 

r A ^bath = Y £ ls(fc|^c|A)s| 2 <5(tofc-^), (47) 

fcEbath 

where \\)s — b\\Q) and |fc)s = b\\0). Here He denotes 
the harmonic coupling between local vibrations and the 
heat reservoir. In fact, after some lengthy algebra, one 
can show that Eq. (|47f is equivalent to the following ex- 
pression: 

TA^bath = — — u^ImlT (lu x ) ux, (48) 
tux 

where ux is the normal mode vector of the Ath local 
vibratio n 48 ! 49 . W(uj), which is defined as the retarded 
heat bath self-energy, is the mechanical counterpart to the 
electronic lead self-energy S r (e). The numerical methods 
used to calculate T, r (e)22r2£. can be directly applied to the 
calculation of n r . 

IV. APPLICATIONS 

In this section we present two applications of our 
ah initio inelastic transport approach. As a first ap- 
plication we calculate inelastic transport properties of 




/t(2Wa) 



FIG. 2. (Color online) Band structure of cumulene for a two- 
atom unit cell. (Red dots: direct DFT calculation; black 
solid lines: Wannier interploted bands), (a) p-type Wannier 
function at an atomic site (b) <r-like Wannier function at a 
mid-bond site. 




Left electrode Device region Right electrode 



FIG. 3. (Color online) Cumulene-benzene-cumulene junction. 



a benzene molecule connected to monoatomic carbon 
chains (cumulenes) . For this benchmark system we apply 
both the SCBA and the LOPT for equilibrium and non- 
equilibrium solutions. As a second example we replace 
cumulene with (3, 3) single- wall carbon nanotube (CNT). 
In this molecular junction we focus on the calculation of 
realistic decay rates and non-equilibrium vibration popu- 
lations. Technically, all DFT calculations are performed 
using the Perdew-Zunger local density approximation 50 , 
norm-conserving pseudopotentials 51 , and a plane-wave 
basis with a cutoff of 55 Ry. 



A. Cumulene-Benzene-Cumulene junction 

We first study a benzene molecule connected to a 
monoatomic carbon chain. It is known that cumulene is 
subject to a Peierls distortion: It readily becomes dimer- 
ized and opens an energy band gap favored by a lower 
energy structure. We thus freeze the structure of cumu- 
lene and use it as a metallic electrode. We construct 
two p-orbitals and one cr-like mid-bond Wannier func- 
tions (see FigfJJ) for this electrode. The two p-orbitals, 
p v and p z , are perpendicular to the transport direction, 
which is along the x-axis in our calculation. Figure [5] 
shows the band structure of cumulene. The energy bands 
are obtained either from a direct planewave-based DFT 
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FIG. 4. (Color online) Differential conductance G = dl jdV 
and its derivative dG/dV with equilibrum vibrationl popu- 
lations calculated using LOPT (black solid line) and SCBA 
(red dashed line). At lower bias two differential conductance 
increases are observed. At higher bias, two large conductance 
drops occur. These conductance changes correspond to peaks 
in dG/dV. 



calculation, or by Wannier band interpolation, and are 
in excellent agreement. While the lowest energy band 
originates from er-orbitals, p-orbitals give rise to doubly 
degenerate 7T-bands around the Fermi level, and trans- 
port properties at the Fermi energy are characterized by 
these two 7T-bands. 

Figure [3] shows the supercell geometry used in the 
transport calculations. The benzene molecule alone is 
allowed to vibrate. The device region, containing the vi- 
brating region and part of the cumulene, is taken to be 
large enough to make sure that there is no direct coupling 
between electrodes, and that the electron-vibration cou- 
pling is zero outside the device region. For the benzene 
molecule, p z -type Wannier functions on carbon atoms 
and cr-like Wannier functions on C — C and C — H bonds 
are constructed. 

The differential conductance G = dl/dV and its 
derivative dG/dV are calculated using either the SCBA 
or the LOPT scheme. Temperature is taken to be 
ksT = lmeV. As seen in FigJH these two approxima- 
tions display essentially perfect agreement. Four con- 
ductance changes are observed in the differential conduc- 
tance curve of FigSJ The corresponding inelastic trans- 
port signals due to electron- vibration interactions appear 
as peaks in the d 2 I/dV 2 plot. The peak position on the 
bias axis corresponds to the vibrational energy involved 
in the electron-vibration scattering events. From four 
peaks observed in FigfJ] one might conclude that there 
are four active vibrational modes, but there is a shoulder 
on the right side of the third peak. This may indicate 
that there is a fifth active vibrational mode. To inves- 
tigate which vibrational modes participate in inelastic 
transport, we performed modernise calculation by keep- 
ing only one particular vibrational mode. For clarity the 
elastic contribution is excluded in the modewise calcula- 
tion. These calculations show that there are five major 
peaks in dG/dV (Fig. [S]). The corresponding vibrational 
configurations of these five active modes are also shown 
in FigJS] While two upward peaks are out-of-plane vi- 
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FIG. 5. (Color online) dG/dV in modewise calculation. Five 
active vibrational modes are found. The corresponding vi- 
brational configuration are illustrated. While the first two 
active modes leading to conductance jumps are out-of-plane 
motions, the three conductance- drop modes correspond to in- 
plane vibrations. 



brations, three downward peaks correspond to in-plane 
motions. 

Thus, electron- vibration interactions can lead to both 
differential conductance rises and drops. This simultane- 
ous occurrence can be understood from scattering theory 
and transmission eigenchannel o 54 ' 55 , following Ref.— . As 
seen in Fig|6j let us consider an electron injected from the 
left electrode on a left-incident ith eigenchannel \^ l L } at 
energy e. If there is no scattering with molecular vibra- 
tions, this electron contributes to elastic conductance 



t(e), 



(49) 



Left -incident ith eigenchannel 

i lTl i / sy transmission 

Vibron emission 
Scattering probability \ ^Tuo 

'■ tjJ 

R 

Right-incident j'th eigenchannel 



FIG. 6. Schematic representation of inelastic scattering in 
the presence of electron-vibration interactions. Solid arrows 
indicate transmission eigenchannels. 




8 



(a) G (2e /h) 
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Bias (V) 
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FIG. 7. (Color online)Inelastic transport calculations 
with nonequilibrium vibrational populations (blue solid line: 
LOPT; red dashed line: SCBA; black dot-dashed line: equilib- 
rium case), (a) differential conductance, (b) second derivative 
of the current, and (c) vibration populations. 



where Go = ^- and Tl^ R {e) is the elastic transmission 
probability of \^ l L (e)). Now let us consider that the elec- 
tron on a left-incident ith eigenchannel \^ l L (e)) is scat- 
tered off to a right-incident jth eigenchannel \^ 3 R {e — fiLS)) 
by emitting a vibrational quanta hui. The outgoing state 
of ^-^(e — hio)) can go either to the left electrode or to the 
right one. Since conductance measured in the right lead 
is considered here, the probability to move back to the 
right electrode for ^^(e — fi^)) is given by its reflection 
probability TZ R (e — fruj). If Vi-^j denotes the probability 
that is scattered to ^^(e — fiw)}, then the total 

conductance can be computed as: 



Gel-vib - Go E 1 ~ E ^-+3 7~l^ R 



+ ]>>^7^(( £ )-M 

ij 

Thus, the conductance change is given by 

AG — G el—vib G no el—vib 



(50) 



E G ° 



.(52) 



Let us examine these relations in detail. First, trans- 
mission and reflection probabilities are approximated by 
those at Fermi energy ef (Since vibrational energies 
are generally small, one might expect that transmission 
and reflection probabilities would not change significantly 
over [ep — hio,eF + huj]). Now, note that conductance 
can increase or decrease depending on the relative mag- 
nitude of Vi^j (tZ r (ef) — 7~l^ R {eF)j , and that while 

(lZ R (eF ) — Tl^.r(ef)j does not depend on a vibrational 

2 



configuration, Vi- 



R^el-vib 



which can be 



calculated from Fermi's golden rul o 56 ' 57 , is determined by 
the electron- vibration interaction matrix. 

For the cumulene-benzene-cumulene system, it is found 
that there are two transmission eigenchannels: the major 
transmission channel R ) with T^ R = T R ^ L = 0.741 
and the minor channel R ) with T^ R = T%^ L = 
0.004. p z orbitals on the cumulene wire and the ben- 
zene molecule constitute the major transmission eigen- 
channel \^ l L R ). p y orbitals on the cumulene wire and a 
bonds of the benzene molecule contribute to the minor 
transmission channel \^> 2 L R ). In addition, while R ) is 
symmetric with respect to the za;-plane, \^ 2 L R ) is anti- 
symmetric. In other words, when P zx denotes the re- 
flection operator with respect to zx-plane, -P Z£C | , I'^ R ) = 

m, R ) and P zx \*l <R ) = -\*% R ) hold. 

The first two active vibrational modes A = 5, 11, shown 
in Figj5j which lead to differential conductance jumps, 
show anti-symmetric vibrational motion with respect to 
the zcc-plane. Then, the corresponding electron-vibration 
interactions 'H^^H satisfy 



H 



A=5,ll 



_ p o/A— 5,11 p 

r ZX iLpl r 7-- 



(53) 



el — vib zx el — vib zx ' 

Because of these reflection symmetries, 

(^T-^i) = <*^£ 6 £l*£> = o, (54) 

1(2) 1(2) 

i.e. scattering from ) to \ty # ) is prohibited. Since 

(K 2 R - T£^ R ) = (R} R - T 2 ^ R ) = 0.255, one obtains a 
differential conductance rise 



AG = Go 



-V. 



2 ■ /v ;,' h. -it'' 

A " n (^-r L 2 ^)l >o. 



2-s-l 



(55) 



In contrast, the last three active modes A = 17, 19 and 
25 are symmetric with respect to the zx-plane: 



t,A=17,19.25 _ p t,A=17,19,25 p 



L el-vib "'^''el-iiiii '21' (56) 

Therefore, one has these reflection selection rules: 



A=17,19,25 hT ,2 \ _ /, T ,2 |-t,A=17,19,25 
e l_ vib |W L ) - \y R \rl el _ mb 



H) = 0- (57) 



Since <p*-\ 7 ' 19 ' 25 > r^ 7 ' 19 ' 25 from numerical calcu- 
lation, then one finds the three differential conductance 
drops 



AG = Go 



7 -,A=17,19,25/t-,1 t-1 x 



^=17.19.25,^2 
+ ' 2->2 



Tl 



L^fRj 



) <0, (58) 



n 



where {ll R -T^ R ) = -0.482 and (H% - , L _> R , 
0.992. More generally this multichannel analysis shows 
that differential conductance rises and drops can occur 
at the same time. 

Next, let us take into account the effect of nonequi- 
librium vibrational populations on transport properties. 
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FIG. 9. (Color online) (3,3) CNT - Benzene - (3,3) CNT su- 
percell geometry used in the decay rate calculations. The vi- 
brating region contains a benzene molecule and three relaxed 
surface CNT layers. 



(a) 



FIG. 8. (Color online)Differential conductance for different 
decay rates (black solid line: h^y\ = 0.1 meV; red solid line: 
hjx = 1 meV; green solid line: tvy\ = 10 meV; blue solid line: 
equilibrium case (fryx — > oo)). 



This situation corresponds to the case where the decay 
rate of a molecular vibration to its surrounding is larger 
than the emission rate due to electron-vibration scatter- 
ing. In order to compare equilibrium and nonequilibrium 
vibration cases the decay rate tvy\ =0.1 meV is chosen 
for all vibrational modes; this condition will be relaxed 
in the next example. 

As shown in FigUJ (a) , nonequilibrium effects lead to 
larger slopes in comparison with the equilibrium case. 
Furthermore, the differential conductance change in- 
creases at the threshold bias voltage. These two changes 
appear as (1) the finite dG/dV value between peaks and 
(2) increased peak heights in FigJTJ (b). When the bias 
exceeds the threshold voltage equal to the vibrational 
energy, the vibrational population starts to increase, as 
observed in FigJT] (c). This is because of the increase in 
phase space of conducting electrons that can emit molec- 
ular vibration quanta. Recalling that electron-vibration 
scattering is roughly proportional to N\, increased vi- 
bration populations enhance inelastic transport signals 
in return. Last, we also considered transport calculations 
where we change the decay rate. Figure [5] shows that the 
differential conductance approaches the equilibrium case 
as the decay rate increases. 
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FIG. 10. (Color online) (a) Decay rates for each vibrational 
mode of the (3,3) CNT-Benzene-(3, 3) CNT junction, (b) 
Decay rate vs. localization (see the text for the definition). 
Decay rates are plotted in a logarithmic scale. 



B. (3,3) CNT - Benzene - (3,3) CNT 

In the previous benchmark, the cumulene wire, sub- 
ject to Peierls' instability, was frozen in order to keep 
its metallic character, and the decay rate for molecular 
vibrations was used as a parameter. Here we replace 
the cumulene wire by a metallic (3,3) carbon nanotubc 
(CNT), which is mechanically stable. Using this elec- 
trode we derive a fully ab initio approach to calculate 
non-equilibrium populations under electron- vibration in- 
teractions. 



Carbon-based nanostructures, such as CNT and 
graphene, could become new platforms for future nan- 
otechnology applications due to their excellent elec- 
tronic properties. Recently carbon-based nanojunc- 
tions have been experimentally fabricated: a pure car- 
bon chain connected to graphene^ or CNTs^, and or- 
ganic molecules coupled to CNT electrodes with amide 
linkages^. These experimental achievements have stim- 
ulated theoretical and computational studies on carbon- 
based nanodevice a 61 ' 62 . In particular, a benzene molecule 
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molecular region 
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FIG. 11. (Color online) Molecular region containing a ben- 
zene molecule, anchoring carbon atoms, and hydrogen atoms 
saturating the CNT edge. 



connected to CNT electrodes (which is the system of our 
interest) was suggested as a molecular switch, operated 
by controlling the relative angle between 7r-orbitals of the 
benzene and the 7r-orbital manifold of CNT electrodes 62 . 
Functionality and performances of molecular devices are 
strongly affected by molecular geometries or anchoring 
points to the electrodes, and they may be affected by vi- 
brations induced by conducting electrons. In the worst 
case, local heating may break down the junctions. 

In our work, we choose a vibrating region by defining 
an extended molecule in which a benzene and the out- 
most relaxed CNT layers are included. This extended 
molecule is seamlessly connected to the bulk CNT elec- 
trodes. The vibrating region contains 56 atoms in total, 
and these correspond to 168 vibrational modes. Figure [9] 
illustrates the supercell geometry used in the decay rate 
calculations: It contains the extended molecule and two 
mechanical principal layers 53 for bulk phonons. 

We proceed as follows. First, by allowing only the 
extended molecule to vibrate, the electron-vibration in- 
teraction Tiei-vib, the vibrational spectrum {hu>\}, and 
the corresponding normal modes are calculated. For the 
decay rates the interatomic force constants for the entire 
supercell in Fig. [9] are calculated. From this interatomic 
force constants one can extract the harmonic coupling 
matrix He between the extended molecule and the bulk 
electrodes. In addition, we take a periodic unit cell of 
the bulk (3,3) CNT and calculate its interatomic force 
constants Hb- 

The calculated decay rates are shown in Fig[TU] (a). 
Note that the decay rates are written in units of eV. 
While most of the decay rates are of the order of 10~ 2 to 
10~ 3 cV, there are few modes with much smaller rates. 
These small decay rates can arise for two reasons. First, 
as seen in FigfTU] (a), decay rates start to increase sig- 
nificantly from the 153th mode on. These vibrational 
modes (between 153 and 168) have energies higher than 
the highest (3,3) CNT phonon energy. Recalling that 
decay processes based on the harmonic coupling essen- 
tially correspond to one vibration to one phonon transi- 
tions, there are no bulk phonons to which the vibrational 
modes lying outside the band width of bulk phonons can 
transfer their vibrational energies. Once the anharmonic 
coupling that makes one vibration to multi- phonons tran- 
sition possible is taken into account, these modes may 
have larger decay rates. However, this correction is be- 
yond our work. 



1: 8.733 meV 
52: 80.48 meV 
77: 102.0 meV 
1 10: 142.4 meV 
149: 188.3 meV 
156: 200.9 meV 



1 0.2 
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FIG. 12. (Color online) Nonequilibrium vibrational popula- 
tions for the most excitable modes as a function of a bias 
voltage. 



Second, notice that among the vibrational modes 
whose energies lie inside the bulk phonon dispersions, 
some of them still have small decay rates. Most of them 
correspond to vibrations that are localized inside the ben- 
zene molecule, or to wagging motions of the surface hy- 
drogen atoms. Because these motions are spatially well 
separated from bulk phonons, one may expect that they 
are less coupled to bulk phonons. In order to measure 
how localized these modes are inside the molecule, we de- 
fine the benzene molecule, two anchoring carbon atoms, 
and the surface hydrogen atoms as the molecular region, 
as seen in FigfTTJ When Vm denotes a projection op- 
erator onto the molecular region, one can find how lo- 
calized the vibration is inside the molecular region from 
|"Pm|A)| 2 where |A) indicates the vibrational state for the 
normal mode A. We call I'PmIA)! 2 the localization mea- 
sure. Figure fT0l (b) shows a relation between the localiza- 
tion measure and decay rates. When the vibration is lo- 
calized in the molecular region, or equivalently I^mIA)! 2 
approaches 1, its decay rate becomes smaller. 

Except for modes between 153 and 168, the majority 
of the vibrational modes overlap with the phonon dis- 
persions of the CNT electrodes. This happens thanks to 
the same chemical character between hydrocarbons and 
carbon-based electrodes. If we were to consider organic 
molecules attached to a metal electrode such as gold or 
platinum, the large mass difference between atoms in the 
molecule and those in the electrode would make most of 
the molecular modes lie outside the electrode phonon dis- 
persions. Therefore, most of the vibrational modes would 
have very small decay rates, significantly increasing the 
probability of the molecular junction to break down. 

Together with Eqs. ([H]), (gU) and (O, we can calcu- 
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FIG. 13. (Color online) Vibrational configurations for the 
modes in Fig |12l 
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late nonequilibrium vibrational populations. Vibrational 
occupations will start to increase as the bias voltage ex- 
ceeds their corresponding threshold voltages. While in 
the vicinity of the threshold voltage the vibrational pop- 
ulations increase linearly, nonlinear effects can appear at 
higher bias voltages. For low-energy modes vibrational 
populations monotonically increase with the bias; how- 
ever, non-monotomic populations are observed for some 
of the high-energy modes. These trends are shown in 
FigH2l where vibrational populations for some of the 
highly excitable modes are illustrated. Mode 1, which 
is low-energy, shows a monotonically increasing behav- 
ior. For the other high-energy modes, their populations 
increase, then decrease in a certain bias voltage range, 
and then start to increase again. 

One can hint at a qualitative explanation for this cool- 
ing behavior examining the local density of states in the 
device region, as recently reported in Refj££. As observed 
in FigfT4l-(a). there is one resonant peak located at a en- 
ergy higher than fi eqi which is the common Fermi level 
before the bias is applied. Figure[T2r(b) shows absorption 
and emission processes where electrons can exploit the 
resonant density of states for different bias voltages. Red 
and blue lines indicate absorption and emission processes 
respectively. For a very small bias none of the electron 
can access the resonant peak, as seen in Figfl4l(b)-(1). 
When the bias increases, the absorption process using 
the resonant peak starts to take place, but the electrons 
participating in the emission process cannot reach the 
resonance (See FigO(b)-(2)). Therefore the absorption 
rate A \ becomes enhanced, so it may lead to a decrease 
in vibrational populations. For higher biases such that 
the resonant peak is located between the left and right 
chemical potentials, the emission process using resonance 





FIG. 14. (Color online) (a) density of states for the device 
region. Close to the equilibrium Fermi level, one resonance 
peak is found, (b) possible absorption (red arrow line) and 
emission (blue arrow line) processes via the resonant peak as 
the bias voltage increases. 



is activated, leading to enhanced emission rates, as shown 
in FigfT4l-(b)-f3). When the bias increases more, another 
resonant emission process becomes possible, as shown in 
Fig03r(b)-(4), while the resonant absorption process at 
which the electrons are reflected back to the left lead is 
prohibited due to Pauli blocking. As a result, the emis- 
sion rate becomes more enhanced in comparison to the 
absorption one, and vibrational populations may increase 
again in this bias range. As an illustratative example, 
Fig ]15l (b) shows the absorption and emission rates for 
mode 156. One can clearly observe that the bias voltages 
for which the absorption and emission rates get enhanced 
are different: The slope of the absorption rate curve first 
increases around 0.4 V, but then decreases at 0.7 V. By 
contrast, the emission rate linearly increases up to 0.8 V, 
and at a bias larger than 0.8 V its slope also increases. 
Difference between bias voltages at which absorption and 
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FIG. 15. (Color online) the absorption (red dashed line) and 
emission (blue solid line) rates for (a) the vibrational mode 1 
(low-energy mode) and (b) 156 (high-energy mode) 
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FIG. 16. (Color online) Total vibrational energy stored in the 
vibrating region, as a function of the mass of the electrode 
atom (black solid line: carbon; red dashed line: silicon; blue 
dot-dahsed line: germanium). 



emission rates become enhanced results in the observed 
cooling behavior in the intermediate bias range. 

For low-energy modes whose energies are much smaller 
than the broadening of the resonant peak, the interme- 
diate cooling regime may not appear distinctly. For ex- 
ample, see in Fig. 1151 (a) the absorption and emission 
rates for mode 1. Unlike mode 156, the absorption and 
emission rates show quite similar dependence on the bias 
voltage, implying that a cooling behavior is not observed. 

Last, not every high-energy mode goes necessarily 
through a cooling down phase, since vibrational popu- 
lations are determined by the interplay of absorption, 
emission, and decay rates. When the decay rate 7a is 
larger than the difference between the absorption and 
emission rates A\ — E\ , the steady state solution can be 



approximated as 

n B (Sw A ) 



7A + A x - E x 



Ex 

7A ' 



(59) 



In this case, the dependence of the vibrational popula- 
tion on the bias voltage becomes similar to that of the 
emission rate, and the cooling behavior does not appear. 

Finally, we would like to stress the importance of mass 
ratios between the conducting molecule and electrodes. 
As pointed out above, since the band width of bulk 
phonons in the electrodes gets narrower as the atomic 
mass of electrodes increases, a molecular junction con- 
nected to electrodes consisting of heavier atoms will have 
less opportunities to thermalize. To demonstrate this 
mass ratio effect, we calculate the total vibrational en- 
ergy stored in the vibrating region by increasing the mass 
of the atoms in the carbon electrode to be e.g. silicon 
or germanium. As shown in Fig |16[ the junction with a 
larger mass ratio has more vibrational energy, and higher 
probability to break down due to heating effects. 



V. SUMMARY 

In this work we have described an ab initio quan- 
tum transport approach based on maximally localized 
Wannier functions to calculate inelastic transport prop- 
erties in the presence of electron- vibration interactions. 
In our implementation calculations are performed using 
the plane-wave DFT code Quantum-ESPRESSO, and 
WANNIER90 for the transformation into MLWF and the 
construction of the Hamiltonian in the Wannier basis. 
Inelastic transport properties such as differential con- 
ductance and non-equilibrium vibrational occupations 
have been calculated within the Meir-Wingreen transport 
formalism based on NEGF. In particular, the electron- 
vibration self-energy is computed using either the SCBA 
or the LOPT. We have tested our implementation ap- 
plying it to carbon-based molecular junctions: a benzene 
molecule connected to cumulene chains and a (3, 3) CNT- 
bcnzene-(3, 3) CNT junction. 

In the first system we have found the differential 
conductance change at bias voltages corresponding to 
the vibrational energies of active modes. Using multi- 
eigenchannel analysis based on scattering theory, we have 
identified that out-of-plane and in-plane vibrational mo- 
tions lead to conductance jumps and drops respectively. 
This analysis on the simultaneous occurrence of conduc- 
tance jumps and drops is consistent with recently re- 
ported calculation^. 

In the second application, where cumulene wires are 
replaced by realistic CNT electrodes, we have focused 
on the decay rates of molecular vibrations and non- 
equilibrium vibrational populations. Small decay rates 
can be rationalized examining the band width of bulk 
phonons that molecular vibration can resonantly access, 
and the localization of molecular vibrations. Also, we 
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have argued that a larger mass ratio between the con- 
ducting molecule and electrodes can enhance local heat- 
ing and ultimately lead to the junction break down. Thus 
organic molecular junctions connected to carbon-based 
electrodes or lighter metals may be more stable in com- 
parison to heavy-metal electrodes (of course, strength of 
the anchoring chemical bonds will also play a role). We 
have also observed that nonequilibrium vibrational oc- 
cupations for high-energy modes can be cooled down in 
an intermediate bias regime, which may be understood 
in terms of a resonant state in the device region. This 
observation agrees well with resonant cooling reported in 
other first-principle calculations 63 . 



These applications and detailed analysis, in a very 
good agreement with other first-principle studies per- 
formed with a local orbital basis set, confirm that our 
ab initio approach to inelastic transport with plane- wave 
basis and Wannier functions could be applied to other 
realistic nanoscale junctions. 
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